#####看727个位点的基因能否富集在DEG数据库中，用基因名富集

##结果：727个位点，共822个基因，300个在DEG中找到，在nobias_AShM中OR=1.19，pvalue=0.02，其他的都不显著

#setwd("E:/0 公共数据库差异情况/db_for_5hmc/CMC等DEG数据库的分析")
setwd("E:/0 公共数据库差异情况/")
library(openxlsx)
sci=read.xlsx("2018 Science supp/aat8127_Table_S1a.xlsx",sheet="DGE")
sci=sci[sci$ASD.fdr<0.1|sci$SCZ.fdr<0.1|sci$BD.fdr<0.1,]
sci=data.frame(symbol=sci$gene_name,ASD.fdr=sci$ASD.fdr,BD.fdr=sci$BD.fdr,SCZ.fdr=sci$SCZ.fdr)

nosva=read.xlsx("CMC/DEG_nosva.xlsx",sheet=1)
nosva=nosva[nosva$adj.P.Val<0.1,]
nosva=data.frame(symbol=nosva$MAPPED_genes,nosva.fdr=nosva$adj.P.Val)

sva=read.xlsx("CMC/DEG_sva.xlsx",sheet=1)
sva=sva[sva$adj.P.Val<0.1,]
sva=data.frame(symbol=sva$MAPPED_genes,sva.fdr=sva$adj.P.Val)

tp_esm=read.xlsx("2018 TP 41398_2018_216_MOESM8_ESM.xlsx",sheet = 2)
tp_esm=data.frame(symbol=tp_esm$X2,tp_esm.fdr=tp_esm$X9)
tp_esm=tp_esm[-1,]
tp_esm$tp_esm.fdr=as.numeric(tp_esm$tp_esm.fdr)
tp_esm=tp_esm[tp_esm$tp_esm.fdr<0.1,]

sig_deg=merge(nosva,sva,by="symbol",all=T)
sig_deg=merge(sig_deg,tp_esm,by="symbol",all = T)
sig_deg=sig_deg[!sig_deg$symbol==".",]		###sig_deg合并完成
length(unique(sig_deg$symbol))
sig_deg=sig_deg[!is.na(sig_deg$symbol),]
sig_deg=sig_deg[-(1:11),]
sig_degid=unique(sig_deg$symbol)

write.csv(sig_deg,"sig_deg.csv",quote=F,row.names = F)

sig_deg=sig_deg[sig_deg$nosva.fdr<0.05|sig_deg$sva.fdr<0.05|sig_deg$tp_esm.fdr<0.05,]
case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
case=tidyr::separate_rows(case,Gene.refGene,sep=";")
caseid=unique(case$Gene.refGene)
caseid=caseid[!caseid=="NONE"]

a=length(intersect(sig_degid,caseid))
b=length(caseid)-a
case1=merge(case,sig_deg,by.x = "Gene.refGene",by.y="symbol")


con.file=list.files(path = "E:/5hmc_file/组织特异性表达/",pattern=".csv")
con.file=con.file[-c(grep("tissue",con.file))]
con.file=paste0("E:/5hmc_file/组织特异性表达/",con.file)

col_names=c("con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=7))
names(result)=col_names
for(i in 1:4){
con=read.table(con.file[i],head=T,sep=",")
con=tidyr::separate_rows(con,Gene.refGene,sep=";")
con$unitID=paste(con$Chr,con$Start,sep = ":")
conid=unique((con$Gene.refGene))
conid=conid[!conid=="NONE"]
c=length(intersect(sig_degid,conid))
d=length(unique(conid))-c
rt_tmp=data.frame(matrix(NA,1,ncol=7))
names(rt_tmp)=col_names
rt_tmp[,1]=con.file[i]
rt_tmp[,2]=a
rt_tmp[,3]=b
rt_tmp[,4]=c
rt_tmp[,5]=d
rt_tmp[,6]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,rt_tmp)
}
 write.csv(result[-1,],"DEG_enrichment.csv",quote=F,row.names = F)